#include "plant_1DSOS_RDM_v3.h"
// --------------------------------------------------
plant_1DSOS_RDM::plant_1DSOS_RDM(double idt,int iDim,int iL,
				 int iMode, double iE_D,
				 double iE_S,double iE_N)
{
  cout << "plant_1DSOS_RDM::plant_1DSOS_RDM()\n";
  // Options --------------------------------------------------------------
  Dim = iDim;
  if(Dim == 1){
    Lx = iL;
    N_CL = 3;
  }else{
    cout << "Dim must be 1";
    exit(-1);
  }

  mode = iMode;
  assert(2*iMode<=iL);
  E_D = iE_D;
  E_S = iE_S;
  E_N = iE_N;

  // 2 Inputs: #0 = T, #1 = W
  int i_num = 2;
  int i_dim[2] = {1,1};
  // 4 outputs:
  //  Output #0: h
  //  Output #1: r2
  //  Output #2: m2
  //  Output #3: a
  int output_num = 4;
  int output_dim[4];
  output_dim[0] = 1;
  output_dim[1] = 1;
  output_dim[2] = 1;
  output_dim[3] = 1;
  block::initialize(i_num,i_dim,output_num,output_dim,idt);

  // States ------------------------------------
  height = new int[Lx];
  neigh = new int[Lx];
  prob_migr = new MCArray(Lx);
  prob_desorp = new MCArray(Lx);

  // Initialize random number generator --------
  gsl_rng_env_setup();
  gsl_T = gsl_rng_default;
  gsl_r = gsl_rng_alloc(gsl_T);
  //  gsl_rng_set(gsl_r,time_seed());
  gsl_rng_set(gsl_r,0);
  double SQRT_TWO = sqrt(2.0);
  double SQRT_PI = sqrt(pi);
  reset();
  // Save a log --------------------------------
  ofstream flog("plant_1DSOS_RDM_v3.log");
  flog << "Lx      = " << Lx << "\n";
  flog << "E_D     = " << E_D << "\n";
  flog << "E_S     = " << E_S << "\n";
  flog << "E_N     = " << E_N << "\n";
  flog.close();
  cout << "END OF plant_1DSOS_RDM::plant_1DSOS_RDM()\n";
}
// --------------------------------------------------
plant_1DSOS_RDM::~plant_1DSOS_RDM(){
  cout << "Begin of ~plant_1DSOS_RDM()\n";
  delete []height;
  delete []neigh;
  delete prob_migr;
  delete prob_desorp;
  gsl_rng_free(gsl_r);
  cout << "End of ~plant_1DSOS_RDM()\n";
}
// --------------------------------------------------
void plant_1DSOS_RDM::reset(){
  //  cout << "Begin of plant_1DSOS_RDM::reset()\n";
  block::time = 0.0;

  // Inputs
  set_temperature(0.0);
  set_ads_rate(0.0);

  // States
  int x;
  // Reset height
  for(x=0;x<Lx;x++){
    height[x] = 0;
  }
  // Reset neighbor count
  for(x=0;x<Lx;x++){
    init_local_stat(x);
  }
  // Update probability
  double site_prob_migr;
  double site_prob_desorp;
  int site_HighNeigh;
  for(x=0;x<Lx;x++){
    site_HighNeigh = neigh[x];
    site_prob_migr = MIGR_RT[site_HighNeigh];
    prob_migr->SetElement(x,site_prob_migr);
    site_prob_desorp = DESORP_RT[site_HighNeigh];
    prob_desorp->SetElement(x,site_prob_desorp);
    migr_ttl = prob_migr->GetTotal();
    desorp_ttl = prob_desorp->GetTotal();
  }
  
  // Outputs
  h  = 0.0;
  r2 = 0.0;
  m2 = 0.0;
  a  = Lx;
  //  cout << "END of plant_1DSOS_RDM::reset()\n";
}
// ---------------------------------------------------
// Count the number of neighbors with larger height
void plant_1DSOS_RDM::init_local_stat(int x){
  neigh[x]=(height[(x+Lx-1)%Lx]>=height[x])+(height[(x+Lx+1)%Lx]>=height[x]);
  return;
}
// ---------------------------------------------------
// Update the table of migration rate table for different number of neighbors
void plant_1DSOS_RDM::init_migr_rate(){
  int i;
  double k_m;
  
  k_m=2*T*KB2H;
  
  for(i=0;i<4;i++){
    MIGR_RT[i]=k_m*exp(-(E_S+i*E_N)*EV2KB/(T+0.0));
  }
  MIGR_RT[4]=0;
  return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_desorp_rate(){
  int i;
  double k_m;
  
  k_m=2*T*KB2H;
  
  for(i=0;i<5;i++){
    //DESORP_RT[i]=k_m*exp(-(E_D+i*E_N)*EV2KB/(temperature+0.0));
    DESORP_RT[i]=0;
  }
  return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::update(double sysTime,double iT,double iW){
  if(time>sysTime)
    return;
	
  double rdm;
  while(time<=sysTime){
    if(iT!=T)
      set_temperature(iT);
    if(iW!=W)
      set_ads_rate(iW);
		
    rate_ttl = ads_ttl+migr_ttl+desorp_ttl;
    rdm = rate_ttl*gsl_rng_uniform(gsl_r);
    if(rdm<ads_ttl)
      adsorption();
    else if(rdm<=migr_ttl+ads_ttl)
      migration();
    else
      desorption();
    
    time += OneMCTime(rate_ttl);
  }
  outputs[0][0] = getH();
  outputs[1][0] = getR2();
  outputs[2][0] = getM2();
  outputs[3][0] = getAutoCorLength();
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::update(double sysTime){
  update(sysTime,inputs[0][0],inputs[1][0]);
}
// --------------------------------------------------
void plant_1DSOS_RDM::getSnapshot(const string FileName){
  ofstream fout(FileName.c_str());
  fout << "Lx   = " << Lx << "\n";
  fout << "Time = " << time << "\n";
  int x;
  for(x=0;x<Lx;x++){
    fout << setw(5) << height[x];
  }
  fout.close();
}
// --------------------------------------------------
double plant_1DSOS_RDM::getH(){
  int x;
  double h_avg,sum = 0;
  for(x=0;x<Lx;x++){
    sum+=height[x];
  }
  h_avg = sum/Lx;
  return h_avg;
}
// --------------------------------------------------
double plant_1DSOS_RDM::getM2(){
  int x;
  double sum = 0;
  double m2;
  double dh;
  for(x=0;x<Lx-1;x++){
    dh = height[x]-height[x+1];
    sum+=dh*dh;
  }
  dh = height[Lx-1]-height[0];
  sum+=dh*dh;
  m2 = sum/Lx;
  return m2;
}
// --------------------------------------------------
double plant_1DSOS_RDM::getR2(){
  int x;
  double R2;
  double sum=0;
  double h_avg;
	
  h_avg = getH();
  for(x=0;x<Lx;x++){
    sum+=(height[x]-h_avg)*(height[x]-h_avg);
  }
  R2 = sum/Lx;
  return R2;
}
// --------------------------------------------------
// Find auto-correleation length. Auto-correlation function is defined
//  in the same way as in Matlab.
//  xcov(m) = \sum_{i=0}^{i=L-m-1}((x(i)-mean(x))*(x(i+m)-mean(x))
// SEE ALSO: Matlab function xcov()
double plant_1DSOS_RDM::getAutoCorLength(){
  int i,m;
  double h_avg = getH();
  double xcov_val_max=0.0;
  // The special case, m = 0.
  m = 0;
  for(i=0;i<Lx-1-m;i++){
    xcov_val_max += (height[i]-h_avg)*(height[i+m]-h_avg);
  }
  // General case, m != 0
  double xcov_val;
  for(m=1;m<Lx-1;m++){
    xcov_val = 0.0;
    for(i=0;i<Lx-1-m;i++){
      xcov_val+=(height[i]-h_avg)*(height[i+m]-h_avg);
    }
    if(xcov_val<0.3679*xcov_val_max){
      break;
    }
  }
  return m*0.2e-9;
}
// --------------------------------------------------
void plant_1DSOS_RDM::set_temperature(double T_set){
  T = T_set;
  init_migr_rate();
  init_desorp_rate();
  k_m2d = exp(-(E_D-E_S)*EV2KB/(T+0.0));
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::set_ads_rate(double W_set){
  W=W_set;
  ads_ttl=W*Lx;
  return;
}
// --------------------------------------------------
// Add or remove an atom from the lattice
void plant_1DSOS_RDM::add_remove_atom(int x, int num){
  int neigh_list[3];
  int neigh_X; 
  int neigh_HighNeigh;
  double site_rate_mig, site_rate_desorp;
  int i;

  // Update height
  if((height[x]==0) && (num==-1)){
    // DO NOTHING if height == 0 and we want to remove an atom
  }
  else{
    height[x]+=num;
  }
  
  GetNeighbors(x,neigh_list);
  // Update neigh counts
  for(i=0;i<N_CL;i++){
    neigh_X = neigh_list[i];
    init_local_stat(neigh_X);
  }
  // Update probability
  for(i=0;i<N_CL;i++){
    neigh_X = neigh_list[i];
    neigh_HighNeigh=neigh[neigh_X];
    site_rate_mig = MIGR_RT[neigh_HighNeigh];
    site_rate_desorp = DESORP_RT[neigh_HighNeigh];
    prob_migr->SetElement(neigh_X,site_rate_mig);
    prob_desorp->SetElement(neigh_X,site_rate_desorp);
    migr_ttl = prob_migr->GetTotal();
    desorp_ttl = prob_desorp->GetTotal();
  }
  return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::GetNeighbors(int x, int *xpos){
  xpos[0] = x;
  xpos[1] = x-1;
  xpos[2] = x+1;
  if(xpos[1]<0) xpos[1]+=Lx;
  if(xpos[2]>=Lx) xpos[2]-=Lx;
  return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::adsorption(){
  add_remove_atom((int)(gsl_rng_uniform(gsl_r)*Lx),1);
  return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::migration(){
  int x;
  int xpos[N_CL];
  int i,n;
  int st[2];
  migr_site(&x);
  if(neigh[x]==4)
    return;
  
  GetNeighbors(x,xpos);
  // Count the number of neighbors with lower height
  for(i=1,n=0;i<N_CL;i++){
    if((height[x]-height[xpos[i]])>0){
      st[n]=xpos[i];
      n++;
    }
  }

  if(n==0) return; // All neighbors are of the same height
	
  // n = (int)(RandN1()*n);
  n = (int)(gsl_rng_uniform(gsl_r)*n);
  add_remove_atom(x,-1);
  add_remove_atom(st[n],+1);

  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::desorption(){
  int x;
  desorp_site(&x);
  add_remove_atom(x,-1);
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::migr_site(int *migr_x){
  double rdm;
  rdm = gsl_rng_uniform(gsl_r)*migr_ttl;
  *migr_x = prob_migr->FindElement(rdm);
  return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::desorp_site(int *desorp_x){
  double rdm;
  rdm = gsl_rng_uniform(gsl_r)*desorp_ttl;
  *desorp_x = prob_desorp->FindElement(rdm);
  return;
}
// --------------------------------------------------
double plant_1DSOS_RDM::OneMCTime(double r_d){
  double rd;
  rd = gsl_rng_uniform(gsl_r);
  while ((rd <= 0) || (rd >= 1)){
    rd = gsl_rng_uniform(gsl_r);
  }
  return -log(rd)/r_d;
}
// -------------------------------------------------
void plant_1DSOS_RDM::getXCov(const string filename){
  int i;
  int m;
  double h_avg = getH();
  double *xcov_val = new double[Lx];
  for(m=0;m<Lx-1;m++){
    xcov_val[m]=0.0;
    for(i=0;i<Lx-1-m;i++){
      xcov_val[m] += (height[i]-h_avg)*(height[i+m]-h_avg);
    }
  }
  
  ofstream flog(filename.c_str());
  for(m=0;m<Lx-1;m++){
    flog << xcov_val[m];
    if(m!=Lx-1){flog << "\n";}
  }
  delete []xcov_val;
  return;
}
